#include "FourierTron.h"
#include <iostream>


using namespace std;
using namespace SeamStress;


std::vector<std::vector<cxfloat> > FourierTron::twiddle;


void FourierTron::calculateTwiddle(unsigned long int N)
{
  unsigned long int NN = N;
  unsigned long int logN = 0;
  while(NN!=0)
  {
    NN = (NN>>1);
    logN++;
  }
  if(twiddle.size()>=logN){return;}
  
  prec <<= 1;
  clearCache(NULL);
  setCachePrec(NULL);
  pins.sewSoftly(&FourierTron::clearCache, nthreads);
  pins.sewSoftly(&FourierTron::setCachePrec, nthreads);
  
  FourierTron fourier(pins.seamstresses, nthreads, prec, cache_size);
  fourier.calculateTwiddleBase(N);
  prec >>= 1;
  clearCache(NULL);
  setCachePrec(NULL);
  pins.sewSoftly(&FourierTron::clearCache, nthreads);
  pins.sewSoftly(&FourierTron::setCachePrec, nthreads);
  
  for(unsigned long int i=0;i<twiddle.size();i++)
  {
    for(unsigned long int j=0;j<twiddle[i].size();j++)
    {
      twiddle[i][j].set_prec(prec);
    }
  }
}


void FourierTron::clearCache(void* arg)
{
  clear_cxfloat_cache();
  clear_refloat_cache();
}


void FourierTron::setCachePrec(void* arg)
{
  init_cxfloat_cache(prec);
  init_refloat_cache(prec);
}


void FourierTron::calculateTwiddleBase(unsigned long int N)
{
  unsigned long int NN = N;
  unsigned long int logN = 0;
  while(NN!=0)
  {
    NN = (NN>>1);
    logN++;
  }
  if(twiddle.size()>=logN){return;}
  
  unsigned long int tsize;
  unsigned long int tsize2;
  
  cxfloat a(0,0,prec);
  
  unsigned long int nsize = 1+nthreads;
  for(unsigned long int i=tempvec1.size();i<nsize;i++){tempvec1.push_back(a);}
  tempvec1[0].set(0,0);
  mul_power_2(pi, 1, tempvec1[0]);
  tempvec1[0].mul_by_i();
  tempvec1[0].neg();
  //now tempvec1[0] = -2*pi*i
  
  
  nsize = twiddle.size();
  
  tempvec2.clear();
  for(unsigned long int i=twiddle.size();i<logN;i++)
  {
    if(i<2){tsize=1;}
    else{tsize = (1<<(i-2));}
    
    for(unsigned long int j=tempvec2.size();j<tsize;j++){tempvec2.push_back(a);}
    twiddle.push_back(tempvec2);
  }
  
  tsize2 = logN-1;
  cxfloat *tbase = &(twiddle[tsize2][0]);
  tsize = (1<<(tsize2+1));
  tsize2 = tsize<<1;
  cxfloat *base = &tempvec1[1];
  (*tbase).set(tsize2,0);
  div(tempvec1[0], (*tbase), (*base));
  exp(*base,*tbase);
  
  for(long int i=(logN-2);i>=((long int)nsize);i--)
  {
    twiddle[i][0].set(twiddle[i+1][0]);
    twiddle[i][0].sqr();
  }
  
  if((twiddle.size() - nsize)<nthreads){thread_tot=(twiddle.size() - nsize);}
  else{thread_tot=nthreads;}
  thread_var1=nsize;
  pins.sewStraight(&FourierTron::twiddleThread1, thread_tot);
  unsigned long int tsize3=(twiddle.size() - nsize);
  twid=0;
  if(nsize==0){twid=1;}
  for(unsigned long int i=twid;i<tsize3;i++)
  {
    twid = (i+nsize);
    tsize = twiddle[twid].size();
    if((tsize-1)<nthreads){thread_tot=(tsize-1);}
    else{thread_tot=nthreads;}
    pins.sewStraight(&FourierTron::twiddleThread2, thread_tot);
  }
  
}


//val is set to exp(-I*j*pi/(2^i))
void FourierTron::twiddleVal(unsigned long int i, unsigned long int j, cxfloat &val)
{
  if(j >= (unsigned long int)(1<<i))
  {
    twiddleVal(i, j - (1<<i), val);
    val.neg();
    return;
  }
  
  while(i>2)
  {
    if((j%2)==1)
    {
      unsigned long int k = ((j-1)>>1);
      unsigned long int s = (1<<(i-1));
      if(k<(s>>2))
      {
        val.set(twiddle[i-1][k]);
      }
      else if(k<(s>>1))
      {
        unsigned long int kk = ((s>>1) - k - 1);
        mpfr_set(val.real_part, twiddle[i-1][kk].imag_part, GMP_RNDN);
        mpfr_set(val.imag_part, twiddle[i-1][kk].real_part, GMP_RNDN);
        val.neg();
      }
      else if(k<((s>>1) + (s>>2)))
      {
        unsigned long int kk = ((s>>1) + k) - s;
        mpfr_set(val.real_part, twiddle[i-1][kk].imag_part, GMP_RNDN);
        mpfr_set(val.imag_part, twiddle[i-1][kk].real_part, GMP_RNDN);
        mpfr_neg(val.imag_part, val.imag_part, GMP_RNDN);
      }
      else
      {
        unsigned long int kk = s-k-1;
        val.set(twiddle[i-1][kk]);
        mpfr_neg(val.real_part, val.real_part, GMP_RNDN);
      }
      return;
    }
    else{j=(j>>1);i--;}
  }
  if(i==2)
  {
    if(j==0)
    {
      val.set(1,0);
    }
    else if(j==1)
    {
      val.set(twiddle[1][0]);
    }
    else if(j==2)
    {
      val.set(0,-1);
    }
    else
    {
      val.set(twiddle[1][0]);
      mpfr_neg(val.real_part, val.real_part, GMP_RNDN);
    }
    return;
  }
  else if(i==1)
  {
    if(j==0)
    {
      val.set(1,0);
    }
    else
    {
      val.set(0,-1);
    }
    return;
  }
  val.set(1,0);return;
}


void FourierTron::twiddleThread1(void *arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  for(unsigned long int i=(thread_var1+w);i<twiddle.size();i+=thread_tot)
  {
    unsigned long int n;
    if((twiddle[i].size()-1)<nthreads){n=(twiddle[i].size()-1);}
    else{n=nthreads;}
    for(unsigned long int j=1;j<n;j++)
    {
      twiddle[i][j].set(twiddle[i][j-1]);
      twiddle[i][j].mul(twiddle[i-1][0]);
    }
  }
  
  
}


void FourierTron::twiddleThread2(void *arg)
{
  unsigned long int w = (*((unsigned long int *)arg));
  twiddleVal(twid+1, thread_tot<<1, tempvec1[1+w]);
  unsigned long int tsize = twiddle[twid].size();
  for(unsigned long int j=(thread_tot+w);j<tsize;j+=thread_tot)
  {
    twiddle[twid][j].set(twiddle[twid][j-thread_tot]);
    twiddle[twid][j].mul(tempvec1[1+w]);
  }
  
  
}


void FourierTron::bitReverse(unsigned long int logN)
{
  vector<unsigned long int> r1;
  r1.assign(2,0);
  vector<vector<unsigned long int> > r2;
  
  for(unsigned long int i=reversal.size();i<=logN;i++)
  {
    r2.clear();
    unsigned long int max = (1<<i);
    for(unsigned long int j=0;j<max;j++)
    {
      unsigned long int j1=j;
      unsigned long int j2=0;
      unsigned long int t1=0;
      unsigned long int shft=0;
      while(j1!=0)
      {
        j2 = j2<<1;
        t1 = j1&1;
        j2 = j2^t1;
        shft++;
        j1 = j1>>1;
      }
      if(shft<i){j2 = (j2<<(i-shft));}
      if(j2<j)
      {
        r1[0] = j2;
        r1[1] = j;
        r2.push_back(r1);
      }
    }
    reversal.push_back(r2);
  }
}

